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Abstract.  We  consider  the  enhancement  of  accuracy,  by  means  of  a  simple  post-processing  technique, 
for  finite  element  approximations  to  transient  hyperbolic  equations.  The  post-processing  is  a  convolution 
with  a  kernel  whose  support  has  measure  of  order  one  in  the  case  of  arbitrary  unstructured  meshes;  if  the 
mesh  is  locally  translation  invariant,  the  support  of  the  kernel  is  a  cube  whose  edges  are  of  size  of  the  order 
of  the  mesh  size  only.  For  example,  when  polynomials  of  degree  k  are  used  in  the  discontinuous  Galerkin 
(DG)  method,  and  the  exact  solution  is  globally  smooth,  the  DG  method  is  of  order  k  +  1/2  in  the  L 2  norm, 
whereas  the  post-processed  approximation  is  of  order  2k  +  1;  if  the  exact  solution  is  in  L 2  only,  in  which  case 
no  order  of  convergence  is  available  for  the  DG  method,  the  post-processed  approximation  converges  with 
order  k  +  1/2  in  L'2  over  a  subdomain  on  which  the  exact  solution  is  smooth.  Numerical  results  displaying 
the  sharpness  of  the  estimates  are  presented. 

Key  words,  post-processing,  finite  element  methods,  hyperbolic  problems 

Subject  classification.  Applied  and  Numerical  Mathematics 

1.  Introduction.  In  this  paper,  we  consider  general  finite  element  methods  for  time-dependent  linear 
hyperbolic  systems  of  the  form 

d. 

ut  +  ^2  Aj  u*j  +  Ao  u  =  0|  (x,  t)  £  ld  x  (0,  T] , 

j=i 

u(x,  0)  =  uo(x),  x  £ 

where  {Aj}d=1  are  real,  constant  coefficient  m  x  m  matrices  such  that  5Zj=i  A j£j  h&s  real  eigenvalues  and 
a  complete  set  of  linearly  independent  eigenvectors  for  all  £  E  Rd ,  and  the  function  u  has  range  in  Wn . 
Our  aim  in  this  paper  is  to  show  how  to  exploit  the  inherently  oscillatory  nature  of  numerical  solutions  to 
this  problem  computed  by  means  of  finite  element  methods  to  enhance  the  quality  of  the  approximation. 
This  enhancement  is  achieved  by  post-processing  the  approximate  solution  only  once,  at  the  very  end  of 
the  computation,  at  t  =  T.  The  post-processing  considered  here  is  completely  independent  of  the  partial 
differential  equation  under  consideration  and  can  be  performed  for  entirely  arbitrary  triangulations;  however, 
it  takes  a  particularly  simple  and  computationally  efficient  form  when  the  triangulation  is  locally  translation 
invariant. 
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To  illustrate  the  basic  idea,  let  us  consider  the  following  simple  model  problem: 

ut  +  ux  =  0,  in  (0, 1)  x  (0,T),  u(x,  0)  =  sin(27r;r)  for  x  E  (0, 1), 

subject  to  periodic  boundary  conditions,  and  let  us  compute  an  approximation  U  to  its  solution  u  by  using 
the  discontinuous  Galerkin  (DG)  method  with  piecewise  polynomials  of  degree  one  over  uniform  grids  of 
spacing  h.  We  also  consider  the  post-processed  approximation  U*  =  I\ff2  *U,  where  the  convolution  kernel 
Kff2(x)  =  j^K4’2(x/h)  is  defined  by 

Ki;2{y)  =  -- ^ip{2)(y  -  1)  +  ^l>(2]{y)  -  j^ip{2)(y  +  1), 

where  is  the  B-spline  obtained  by  convolving  the  characteristic  function  =  y  of  the  interval 
(  —  1/2, 1/2)  with  itself  once.  In  Fig.  1.1  we  display,  for  T  =  0.1  and  h  =  1/10  and  h  =  1/20,  the  er¬ 
rors  x  i — }  u(T,x )  —  U(T,x )  and  x  H-  u(T,x )  —  U*(T,x).  The  time-step  was  chosen  so  small  that  the  overall 
accuracy  of  the  method  is  dominated  by  the  spatial  error.  We  note  the  oscillatory  nature  of  the  error 
x  u(T,x )  —  U(T,x )  typical  of  finite  element  methods  and  the  apparent  superconvergence  of  the  numerical 
solution  at  the  two  Gauss-Radau  points,  a  fact  discovered  in  1995  by  Adjerid,  Aiffa,  and  Flaherty  [2];  see 
also  their  recent  work  [1],  In  contrast  with  this  behavior,  we  observe  the  complete  absence  of  oscillations 
from  the  error  u(T)  —  U*(T).  This  shows  that  convolving  the  approximate  solution  U  with  the  kernel  K4f2 
filters  out  the  numerical  oscillations  around  the  exact  solution.  Moreover,  the  result  of  such  a  filtering  is 
a  new  approximation  U*  that  converges  faster  to  u  than  U.  Indeed,  in  Fig.  1.2,  we  display  the  functions 
x  log(  |  u(T,  x)  —  U(T ,  x)  | ),  for  h  =  1/10, 1/20, 1/40  and  1/80;  we  observe  that  each  time  h  is  halved,  the 
maximum  of  x  \u(T,x)  —  U*(T,x)  |  is  divided  by  a  factor  not  less  than  eight.  This  indicates  that  the 
post-processed  approximation  is  at  least  third-order  convergent;  the  original  approximate  solution  U  exhibits 
only  second-order  convergence. 

In  Figs.  1.3  and  1.4  we  repeat  the  above  experiment  using  polynomials  of  degree  two.  Again  we  observe 
the  oscillatory  nature  of  the  approximation  and  the  superconvergence  at  the  three  Gauss-Radau  points  in 
Fig.  1.3  (top),  and  that  the  oscillations  are  filtered  out  upon  convolution  in  Fig.  1.3  (bottom).  This  time, 
the  convolution  kernel  K^’3(x)  =  j^K6’3(x/h)  is  defined  by 

K6,3(y )  —  2)  — >^—ibt'3\y  —  1)  —  ——,ipt'A\y) 

1920  ’  480  ’  320  ’ 

— (y  +  1)  H — (y  +  2) , 

480 V  ’  1920  ’ 

where  ip (4)  is  the  B-spline  obtained  by  convolving  the  characteristic  function  t/i(1)  =  y  of  the  interval 
(  —  1/2, 1/2)  with  itself  three  times.  In  Fig.  1.4,  we  see  that  each  time  h  is  halved,  the  maximum  error 
decreases  by  a  factor  not  less  than  thirty  two.  This  shows  that  the  error  in  the  post-processed  approxima¬ 
tion  is  of  fifth  order. 

In  connection  with  this  fact,  we  note  here  that  in  1996  Lowrie  [16]  found  analytical  and  numerical  evidence 
that  when  polynomials  of  degree  k  are  used,  a  ‘component  of  the  error’  of  the  DG  method  converges  with 
order  2  k  +  1  in  the  L2  norm;  this  fact  stands  in  striking  contrast  with  convergence  of  order  k  +  1/2  for 
the  underlying  DG  approximation  (k  +  1  for  the  one-dimensional  case  and  special  grids  in  several  space 
dimensions).  In  this  paper,  we  provide  a  firm  mathematical  basis  for  this  observation,  and  show  how  to 
compute  the  superconvergent  approximation  U*  by  a  simple  post-processing  technique  which  is  independent 
of  the  equation  and  of  the  numerical  method. 
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Fig.  1.1.  The  errors  u  —  U  (solid  line)  and  u  —  U*  (dots)  at  T  =  0.1  for  h  =  1/10  (top)  and  h  =  1/20  (bottom).  The 
function  u  is  the  smooth  exact  solution,  U  is  the  approximation  given  by  the  DG  method  with  polynomials  of  degree  one,  and 
U*  =  K't2  *  U. 


Fig.  1.2.  The  errors  log(|  u  —  U*  |)  atT—  0.1  for  h  =  1/10  (top),  h  =  1/20,  h  =  1/40,  and  h  =  1/80  (bottom).  Each  time 
h  is  halved,  the  maximum  error  decreases  by  a  factor  not  less  that  8;  the  order  of  convergence  is,  therefore,  not  less  than  3. 


The  paper  is  organized  as  follows.  In  Section  2,  we  present  a  brief  account  of  the  development  of  the  ideas 
behind  this  paper.  In  Section  3,  we  state  and  discuss  our  main  theoretical  results,  and  in  Section  4  we  present 
their  proofs.  In  Section  5,  we  display  numerical  experiments  which  not  only  verify  our  theoretical  results  but 
also  indicate  how  this  kind  of  post-processing  can  be  applied  to  convection-diffusion  and  non-linear  problems. 
We  conclude,  in  Section  6,  with  some  remarks. 

2.  A  brief  overview  of  the  development  of  post-processing  techniques.  In  order  to  introduce 
the  basic  ideas  of  our  work  and  to  put  them  into  proper  perspective,  we  briefly  review  the  development  of 
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Fig.  1.3.  The  errors  u  —  U  (solid  line)  and  u  —  U*  (dots)  at  T  =  0.1  for  h  =  1/10  (top)  and  h  =  1/20  (bottom).  The 
function  u  is  the  smooth  exact  solution,  U  is  the  approximation  given  by  the  DG  method  with  polynomials  of  degree  two ,  and 


U*  =  Jc  ®’3  *  U. 


Fig.  1.4.  The  errors  log(|  u  —  U*  |)  atT—  0.1  for  h  =  1/10  (top),  h  =  1/20,  h  =  1/40,  and  h  =  1/80  (bottom).  Each  time 
h  is  halved,  the  maximum  error  decreases  by  a  factor  not  less  that  32;  the  order  of  convergence  is,  therefore,  not  less  than  5. 


post-processing  techniques  devised  to  improve  the  quality  of  numerical  approximations.  For  further  details, 
the  reader  should  consult  the  monograph  of  Wahlbin  [21]  on  superconvergence  in  Galerkin  finite  element 
methods. 

2.1.  Finite  difference  and  spectral  methods  for  hyperbolic  problems.  In  1977,  Majda  and 
Osher  [18]  considered  formally  high-order  accurate  dissipative  difference  schemes  for  hyperbolic  problems. 
They  studied  a  one-dimensional  model  problem  of  a  two-by-two  hyperbolic  system  whose  characteristics  are 
parallel  to  x  =  ±t;  the  initial  condition  is  a  step  function  whose  discontinuity  is  located  at  the  origin.  Majda 
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and  Osher  showed  that  the  rate  of  convergence  on  the  region  between  the  characteristics  issuing  from  the 
origin,  x / 1  <  1  —  S'2,  is  independent  of  the  numerical  scheme.  They  pointed  out  that  in  1962  Fedorenko 
[11]  and  in  1969  Apelkrans  [3]  displayed  numerical  evidence  that  the  order  of  convergence  had  to  be  one. 
However,  by  selecting  a  suitable  approximation  of  the  initial  datum,  Majda  and  Osher  showed  that  the  order 
of  convergence  can  be  increased  to  two.  Moreover,  they  found  that  they  could  recover  the  full  formal  order 
of  accuracy  of  the  scheme  on  the  region  \x/t\  <  1  —  S'2  provided  they  preprocessed  the  initial  data  in  an 
appropriate  way.  In  1986,  Johnson  and  Pitkaranta  [14]  used  a  similar  idea  in  the  analysis  of  the  DG  method 
for  linear  hyperbolic  problems.  The  question  of  post-processing  the  initial  data  is  considered  in  the  book 
of  Brenner,  Thomee  and  Wahlbin  [5];  see  also  the  work  of  Jovanovic,  Tvariovic  and  Siili  [15]  concerning  the 
use  of  convolution  mollifiers  with  B-spline  kernels  for  second-order  hyperbolic  boundary  value  problems  with 
non-smooth  data. 

In  1978,  Mock  and  Lax  [19]  showed  that  for  a  difference  scheme  of  any  formal  order  of  accuracy  p.,  for 
linear  hyperbolic  systems,  the  moments  of  the  exact  solution  converge  with  order  p,  provided  that,  again,  the 
initial  data  was  suitably  preprocessed.  This  result  holds  even  if  the  exact  solution  contains  discontinuities. 
They  also  showed  how  to  post-process  the  approximate  solution  by  a  simple  convolution  to  enhance  its 
accuracy  over  regions  of  smoothness  of  the  exact  solution:  if  the  solution  was  sufficiently  smooth  locally, 
they  could  obtain  nearly  the  full  order  of  convergence  p.  provided  that  the  support  of  the  kernel  was  of  order 
almost  one.  This  seems  to  have  been  the  first  instance  when  the  ideas  of  (i)  preprocessing  the  initial  data, 
(ii)  obtaining  error  estimates  for  the  moments,  and  (iii)  post-processing  the  approximation,  appear  clearly 
delineated. 

Later,  in  1985,  Gottlieb  and  Tadmor  [12],  motivated  by  the  work  of  Mock  and  Lax  [19],  found  a  spectrally 
accurate  post-processing  kernel  for  spectral  methods;  see  also  the  1978  paper  by  Majda,  McDonough  and 
Osher  [17].  Again,  the  full  spectral  accuracy  could  be  recovered  by  using  a  convolution;  the  measure  of  the 
support  of  the  kernel  had  to  be  of  order  one. 

2.2.  Finite  element  methods  for  elliptic  problems.  Quite  independently  of  the  developments 
reviewed  above,  in  1977  Bramble  and  Schatz  [4]  considered  linear  elliptic  problems  and  showed  how  to 
post-process  the  finite  element  solution  by  means  of  a  simple  convolution  to  enhance  the  quality  of  the 
approximation.  They  showed  that  the  order  of  convergence  could  be  doubled  if  the  exact  solution  was 
locally  smooth.  It  is  important  to  point  out  that,  just  like  Mock  and  Lax,  Bramble  and  Schatz  proved  a 
negative-order  norm  error  estimate  (an  error  estimate  of  the  moments  in  Mock  and  Lax’s  terminology)  and 
then  showed  how  to  use  it  to  enhance  the  approximation  by  a  convolution.  However,  unlike  Mock  and  Lax’s 
convolution  kernel,  for  locally  translation  invariant  grids  the  Bramble-Schatz  kernel  has  support  in  a  cube 
whose  diameter  is  of  order  h  only;  this  fact  represents  a  considerable  advantage  from  the  computational 
point  of  view. 

Also  in  1977,  Thomee  [20]  extended  the  work  of  Bramble  and  Schatz  [4]  to  include  superconvergence  of 
the  derivatives  and  gave  an  elegant  proof  of  their  approximation  results  by  using  Fourier  analysis. 

An  application  of  the  Bramble  and  Schatz  technique  to  the  simulation  of  miscible  displacement  was 
devised  and  analyzed  by  Douglas  [9];  other  applications  can  be  found  in  the  book  of  Wahlbin  [21]. 

2.3.  The  main  ideas.  In  this  paper,  we  apply  the  ideas  of  Mock  and  Lax  [19]  and  Bramble  and  Schatz 
[4]  to  enhance  the  accuracy  of  finite  element  approximations  to  hyperbolic  problems  by  post-processing. 
We  proceed  as  follows.  First,  we  obtain  an  estimate  of  the  error  between  the  analytical  solution  u  and  the 
post-processed  numerical  approximation  U  in  terms  of  negative-order  Sobolev  norms  of  u  —  U.  This  result 
does  not  depend  on  the  partial  differential  equation  under  consideration  or  on  the  numerical  scheme.  Next, 
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we  obtain  negative-order  norm  a  priori  estimates  for  the  error  between  the  exact  solution  of  a  hyperbolic 
problem  and  its  finite  element  approximation  U.  The  final  error  estimate  is  then  obtained  by  combining  the 
above  bounds. 

3.  The  results.  In  this  section,  we  present  and  discuss  our  main  theoretical  results. 

3.1.  An  approximation  result.  We  begin  by  presenting  a  result  that  relates  negative-order  norm  a 
priori  estimates  of  the  difference  between  u  and  an  arbitrary  approximation  U  for  u  to  L2-error  estimates  of 
the  difference  between  u  and  the  post-processed  counterpart  U. 

Let  us  recall  the  definition  of  a  negative-order  Sobolev  norm  on  an  open  set  0  C  We  denote  by 
II  u  ||o,o  the  standard  T2-norm  of  u  on  f l.  For  any  natural  number  l ,  we  consider  the  norm  and  seminorm  of 
the  Sobolev  space  defined  by 


r 

^  1 1/2 

( 

^  11/2 

u  Hi, a  =  \ 

E  II  D°u  ll5,o  }  * 

u  =  < 

E  II  D°u  llo.o  } 

l 

\a\<C  J 

l 

|a|  =  f  J 

Sobolev  norms  and  seminorms  for  vector-valued  functions  from  17^(0,  Mm)  are  defined  analogously  and  are 
denoted  by  the  same  symbol  as  in  the  scalar  case.  We  then  define  the  negative  order  Sobolev  norm  |j  • 
t  >  1,  by 

fQu(x).cf>{x)dx 

II  u  ll-po  =  sup  ,,  — ,, - • 

<bec^(Q)  110  111, n 

Negative-order  norms  can  be  used  to  detect  the  oscillations  of  a  function  around  zero.  For  example,  for 
fi  =  (  —  1, 1),  l  >  1  and  un(x)  =  sin(27r  N  x),  a  simple  computation  gives  ||  ujv  ||-f,Q  =  1/(2 ir  N)e,  indicating 
that  un  oscillates  about  zero  in  a  very  regular  manner. 

Next,  we  describe  the  type  of  post-processing  to  be  considered  following  Bramble  and  Schatz  [4].  We 
post-process  the  approximate  solution  by  convolving  it  with  a  kernel  K^(x)  =  Kv,e (x / H) / Hd  which  has  to 
satisfy  three  properties;  the  first  of  these  is  that  Kv,e  has  compact  support.  The  second  is  that  it  reproduces 
polynomials  p  of  degree  v  —  1  by  convolution,  that  is, 

Kv,t  *  p  =  p. 

This  is  the  type  of  kernel  used  by  Mock  and  Lax  [19].  The  kernels  used  by  Bramble  and  Schatz  [4]  which 
we  shall  next  describe  have  the  further  property  that  they  are  linear  combinations  of  B-splines.  Let  \  be 
the  characteristic  function  of  the  interval  (  —  1/2, 1/2)  and  let  S  denote  the  Dirac  distribution  concentrated 
at  x  =  0.  Then,  we  define  recursively  the  functions  as  follows: 

^p(O)  =  S _  ^(n+l)  =  jp(n)  *  Xj  for  n  >  0, 

and,  given  an  arbitrary  multi-index  a  =  (op  , . . . ,  ad)  and  y  =  (ylif . .  ,yd)  £  Md,  we  set 

'.'(Ol'0yi )  •  •  •  (//,/)- 

We  also  set  1  =  (1, ... ,  1).  The  third,  and  final,  property  of  the  kernels  considered  here  is  that  they  are  of 
the  form 


Kv’n(v)=  E  7), 

7  eX* 


(3.1) 
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where  k”’n  £  M.  Note  that  since  the  support  of  Kv’n  has  been  assumed  compact,  there  are  only  finitely 
many  non-zero  coefficients  kT’n  in  this  sum. 

The  imposition  of  these  hypotheses  is  motivated  by  the  following  observations:  the  compactness  of  the 
support  of  the  convolution  kernel  is  advantageous  from  the  computational  point  of  view;  the  second  property 
ensures  that  accuracy  of  order  v  is  not  destroyed  by  post-processing;  the  third  property  allows  us  to  express 
derivatives  of  the  convolution  with  the  kernel  in  terms  of  simple  difference  quotients.  Indeed,  it  is  very  easy 
to  verify  that  for  multi-indices  a  and  3  such  that  3i  >  a*  for  |  =  1, . . . ,  d,  we  have 


£>a(V4f  *u)  =i^-a>  *d%v 


—  jSd-w) 


(3.2) 


where  tp^ix)  =  ip(-a/H'1  /Hd, 

dH  ■=  dH,i  ■  ■  ■  dHd,d  and  dHjv(x)  =  ^(v{x  +  ^Hej)-v(x-^HeJ)). 

This  fact  can  then  be  exploited  in  the  finite  element  framework,  as  will  be  seen  later.  We  are  now  ready  to 
state  an  approximation  result  which  shows  that  local  smoothness  of  u  on  the  one  hand  and  negative-order 
norm  estimates  of  divided  differences  of  the  error  u  —  U  on  the  other  lead  to  a  local  bound  on  u  —  K'{jn  *  U 
in  the  L2-norm. 

Theorem  3.1  (Bramble  and  Schatz[4]).  Let  v  and  t  be  two  natural  numbers.  Suppose ,  further ,  that 
Kfjll(x)  =  Kv’n (x / H) / Hd  where  Kv,n  is  a  function  of  compact  support  which  reproduces  polynomials  of 
degree  v  —  1  by  convolution,  and  which  is  the  linear  combination  of  B-splines,  as  in  (3.1).  LetU  be  a  function 
in  L2(fli),  where  Qi  is  an  open  set  in  W1' ,  and  let  u  be  a  function  in  IF(fli).  Let  Slo  be  an  open  set  in 
such  that  fl0  +  2supp(A'^n)  CC  Oi  for  all  H  <  Ho.  Then,  for  H  <  Ho,  we  have 

||  u  -  K;,'1  *  U  ||0,n0  <  ^  Ci  \u  U.q,  +  Ci  C-2  Y,  II 9h  («  -  U)  ||_^,ni , 

M<t 


where  C\  =  Yl~,£Zd  I  \  and  C-2  depends  solely  on  flo,  Oi,  d,  v ,  and  l. 

To  illustrate  the  importance  of  this  result,  let  us  assume  that  there  exist  real  numbers  p.  >  0  and  a  £  [0,  £] 
such  that,  for  all  H  <  H0, 


Y,  \\d%  (u  -  U)  II-eq,  <  C3  h»  H~a.  (3.3) 

\a\<l 

Note  that  the  number  a  measures  how  well  it  is  possible  to  estimate  the  negative-order  norm  of  the  divided 
differences  of  u  —  U.  In  the  worst  case,  a  =  l\  this  is  the  case  treated  by  Mock  and  Lax  [19].  In  the  finite 
element  framework,  however,  it  is  possible  to  take  a  to  be  different  from  l,  as  Bramble  and  Schatz  [4]  showed 
for  second-order  elliptic  problems. 

Inserting  the  inequality  (3.3)  in  the  inequality  of  Theorem  3.1,  we  get 

||  u  -  Kfjn  *  U  Ho.no  <  ^  C-i  |  u  |„,ni  +  Ci  C2  C3  H~a 

<  Ci  max{ |  u  v.uJiA.  Ci  C3 }  {Hv  +  h?  H~a). 


If  we  now  define  H  to  be  the  solution  of  the  equation  Hv  =  /d'  H  a,  we  obtain  the  following  result. 

Corollary  3.2.  Let  the  hypotheses  of  Theorem  3.1  hold,  and  suppose  that  (3.3)  is  valid.  Then,  for 
H  =  /d'/C+C  <  H0,  we  have 

II  u  —  Kv^n  *U  ||0,n0  <Che\ 
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where  C  =  2C\  max{|  u  \v,o.i  / v\,C-2  C 3}  and  6  =  v/(y  +  a). 

Note  that  in  the  worst  possible  case,  that  is  when  a  =  |,  this  implies  that 

\\u-K".n*U \\o,n0 

with  0  =  v/{  v  + 1)  <  1.  The  only  possibility  we  then  have  for  raising  the  order  of  convergence  is  to  hope 
that  the  function  u  is  very  smooth  so  that  we  can  choose  v  large  and  positive.  Unfortunately,  even  if  this 
were  actually  possible,  the  support  of  the  convolution  kernel  would  be  contained  in  a  cube  whose  diameter 
is  of  order  H  =  which  converges  to  a  quantity  of  order  one  as  v  increases  to  infinity;  this  in  turn 

renders  the  evaluation  of  the  convolution  computationally  inefficient. 

On  the  other  hand,  in  the  best  possible  case  (that  is  when  a  =  0),  taking  v  =  p,  would  permit  choosing 
0  =  1,  H  =  h  and  we  would  then  have 

||  u  -  Kv/*  -kU  ||o,o0  <Chr 

In  other  words,  for  a  =  0  we  obtain  the  same  order  of  convergence  for  u  —  K *  U  in  the  local  L2  norm 
as  that  of  the  local  negative-order  norm  error  estimate  in  (3.3).  Moreover,  this  is  achieved  by  using  a 
convolution  kernel  whose  support  is  contained  in  a  cube  whose  diameter  is  of  the  order  H  =  h  only;  this 
renders  the  evaluation  of  the  convolution  a  very  fast  computation.  The  examples  shown  in  the  Introduction 
correspond  to  this  case  with  v  =  p,  =  2k  +  1,  where  k  is  the  degree  of  polynomials  in  the  discontinuous 
Galerkin  method. 

3.2.  Negative-order  norm  error  estimates  for  finite  element  methods. 

3.2.1.  The  weak  solution.  As  stated  in  the  Introduction,  we  consider  the  following  Cauchy  problem: 

d 

ut  +  ^2  +  Aou  =  0,  (x,t)  rfx  (0,  T],  (3.4) 

j= 1 

u(x,  0)  =  uo(x),  x  €  Md,  (3.5) 

To  make  the  presentation  of  the  ideas  as  simple  as  possible,  we  reduce  unessential  technicalities  by 
assuming  that  the  matrices  in  the  equation  (3.4)  are  independent  of  time  and  space  and  by  taking  the  initial 
data  to  be  1-periodic  in  each  of  the  coordinate  directions  Xi,i  =  1, . . . ,  d,  and  we  seek  a  solution  to  the  above 
problem  which  is  1-periodic  in  each  coordinate  direction. 

We  suppose  that  the  system  of  equations  (3.4)  is  strongly  hyperbolic,  that  is,  there  exists  a  family  of 
real  m  x  m  matrices  {£(£)  :  £  €  ld}  and  a  constant  K  >  0  such  that 

j= 1 

is  a  diagonal  matrix  for  all  and 

sup  (||  S(0  ||  +  ||  S-1(£)  || )  <  K.  (3.6) 

l€l=i 

Letting  I  =  (0,  l)d,  the  weak  solution,  u(x,t),  of  (3.4)  satisfies 

i  d 

(u,ip)i(t)  =  (uo,ip(0))i  +  /  (u,ipt  +  ^2  A*jVxj  -  A*0v)i  dr 

0  j  1 
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(3.7) 


for  all  ip  £  C°°  ( [0,  T] ;  Hpel{Rd ,  Km ))  and  t  £  [0,  T]  where  A*  is  the  transpose  of  Aj  and  in  the  above  equation 
and  below 

(u,tp)i(t)  =  j  u(x,t)p>{x,t)  da;. 

Here,  Hpei(Rd ,Wn)  denotes  the  Sobolev  space  of  1-periodic  functions  defined  as  follows.  Let  C“r(Kd, Km)  be 
the  subset  of  C°°  ( Rd ,  Rm )  of  1-periodic  functions .  We  then  define  Hpel  ( Rd ,  Rm )  as  the  closure  of  ( Rd ,  Rm ) 
for  the  H1(I,Rm)~ norm. 

It  follows  from  (3.6)  that  the  problem  (3.7)  is  well  posed  in 

L;er(Kd;  Km)  =  {/  £  L~oc(Rd;Rm)  :  f(x  +  a)  =  f(x)  for  all  x  £  Md,  a  £  Zd} 

with  respect  to  the  norm  ||  •  \\l2(I)',  see  Theorem  6.3.2.  on  p.  219  of  [13]. 

3.2.2.  The  finite  element  methods.  Next,  we  describe  the  class  of  finite  element  approximations 
to  (3.4).  It  includes  the  standard  Galerkin  method,  the  Galerkin  method  with  artificial  diffusion  and  the 
discontinuous  Galerkin  method.  With  slight  modifications  we  could  have  easily  included,  for  example,  the 
streamline  diffusion  method  and  the  stabilized  discontinuous  Galerkin  method;  however,  in  order  to  avoid 
unnecessary  technical  complications,  we  have  chosen  not  to  consider  these. 

Let  7*  =  {  K}  be  a  regular  triangulation  of  Rd ,  invariant  under  translations  by  a  £  Zd,  whose  elements 
I\  are  open  and  have  diameter  Ii  k  less  than  or  equal  to  h.  It  will  be  assumed  throughout  that  each  I\  £  % 
is  contained  either  in  I  or  in  Rd  \  I.  For  a  nonnegative  integer  k,  we  associate  with  the  triangulation  %  the 
broken  Sobolev  space 

tfper,ft(«d; «ro)  =  TlKerhHk(K-, Km)  n  L;ei(Rd- Mm). 

For  k  =  0,  we  shall  write  Lperh(Rd ;Rm)  =  Hpelh(Rd;Rmj,  We  then  consider  two  finite  element  sub¬ 
spaces  Mh  and  A 4.  of  Hpel  h(Rd',Rm),  and  the  broken  Lo  inner  product  ( • ,  • ) ft  defined  on  Lpel  h(Rd  ;Rm)  x 

Ller,h(Rd’Rm)  by 

(W,x)h=  E  (W,x)k,  (3.8) 

KeThiI 

where  Th,i  =  {K  £  T  :  K  C  I}. 

We  define  the  finite  element  approximation  U  :  [0,T]  — y  Mh  as  the  solution  to 

(Ut(t),x)h  +  B(U(t ),  x)  =0,  x  £  A fh,  (3.9) 

U(0)  =  Phuo,  (3.10) 

where  B(-,-)  is  a  bilinear  form  defined  on  Mh  x  iLper  h(Rd;Rm),  and  the  operator  Ph  :  iper(Kd;  lRm)  — >  Mh 
is  the  orthogonal  projection  in  the  norm  of  L'2(I). 

In  Table  3.1,  we  describe  different  choices  of  the  form  B  that  give  rise  to  different  finite  element  methods; 
in  each  of  these  Mh  =  Afh,  although  this  need  not  be  the  case  in  general. 

The  operator  A  and  the  bilinear  form  {•,  •)/,  that  appear  in  Table  3.1  are  defined  as  follows: 

d 

-4.x  =  E  ' \  .£  +  !u\. 

3=  1 

d 

-4*x  =  -  E  AjXgj.  +  ir.v- 

j=i 

(A  U,  x)  ft  =  E  E  n  U  ’  7)  e , 

KeTh.i  eedK 


9 


Table  3.1 

Examples  of  finite  element  methods. 


Method 

M  h  C  C° 

B(U,V) 

Standard  Galerkin  (SG) 

yes 

[AU, if)! 

SG  with  artificial  diffusion 

ves 

(AU,  /?)/  +  h1  (VG,  Vi]) i ,  7>  1, 

Discontinuous  Galerkin 

no 

(U,A*v)h  +  {AU,ri)h 

where  A  ■  n  =  A\ni  +  •  •  •  +  Ajnd,  n  =  (ni5 . . . ,  rid)  is  the  unit  outward  normal  vector  to  K  on  e  C  dK ,  and 
U  is  the  numerical  flux  of  the  DG  method  defined  as  follows.  Given  an  element  K  and  a  face  e  €  d K,  let  us 
denote  by  Ke  G  7/,  the  element  sharing  the  edge  e  with  I\  and  denote  by  Uk  and  Uk.;  the  traces  of  U  on  e 
from  I\  and  Ke,  respectively.  We  compute  the  m  x  m  diagonal  matrix  diag(Ai , . . . ,  Am)  =  S(n)  (A-n)  5_1  (n) 
and  set  V  =  5_1('/r)  U  and 


Vj  = 


( Vk)j 
(VkJj 


The  numerical  flux  is  defined  as  follows: 


if  A j  >  0, 
otherwise  . 


0  =  S(n )  V. 


(3.11) 


3.2.3.  The  negative-order  error  estimate.  We  now  give  sufficient  conditions  for  the  finite  element 
method  which  ensure  that,  for  a  given  time  T,  our  approximate  solution,  U(T),  converges  with  high  order 
in  a  negative-order  norm  over  a  given  subdomain  O0  CC  I  to  the  weak  solution  u(T).  Given  that  l  >  0,  we 
wish  to  estimate 


«(T)-£F(T)||_^b 


sup 

$GC^°(Q0) 


{u(T)-U(T),$) 

II  ^  I  It, 0,o 


We  begin  by  considering  the  solution  to  the  dual  problem:  Find  a  function  such  that  <p(-,t)  is  1-periodic 
in  each  coordinate  direction  for  all  t  G  [0,T)  and 


d 


^  +  E  A*J  'Axi  -  A  =  0, 

j  =  1 

in  Rd  x  (0, T), 

(3.12) 

<p(x,T)  =  $(x), 

x  eRd, 

(3.13) 

where  $  is  an  arbitrary  function  in  C'o°(flo). 


{u(T)  -  U(T),$)  =  (u,  <p)(T)  -  (U,W)(T) 

(«u-v(()))  (r:^)(T) 

=  (uo,<p(0))  -  (Ut<p)(0)  +  J  ^ ( U,  tp)  dr 
=  (uo  ~  PhUo,<p( 0))  -  /  {(Ut,ip)  +  (U,tpt)}  dr. 
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Since,  by  (3.9),  for  x  :  [0,T]  ->•  Afh, 

I  (/',./)  dr  /  (Ut,<p  —  x)  dr  +  /"  (Af,x)dr 
«/o  </o  «/o 

=  /  (t4,y>  -  x)dr  -  [  B(U,x)dr 
Jo  Jo 

=  /  -  x)  +  B(U,<p  -  x)}  dr  - 

Jo 


B(U,V)&r, 


we  obtain  that 


(u(T)  -  U(T),  $)  =  0M  +  0jv  +  0c,  (3.14) 

where 

0  m  =  (u0  -  PhU0,ip( 0)), 

Qn  =  ~  f  {( Ut ,  vJ  V)  +  B(U,  ip  -  x)}  dr, 

Jo 

Qc  =  -[  {(U,<pt)-B(U,<p)}  dr. 

Jo 

Next,  we  introduce  some  general  assumptions  on  Ad/j  and  A4  which  will  enable  us  to  estimate  these  three 
terms. 

Let  Do  CC  Hi  C  I,  r  >  0,  £  >  1,  and  suppose  that  uo  £  Lper(K'<;  Km)  fl  Hr{VPl\ ;  Km),  where  "Dili 
denotes  the  domain  of  dependence  for  the  set  fh;  see  Fig.  3.1. 


Fig.  3.1.  Example  of  the  domain  of  smoothness  of  u(T),  fig,  and  of  a  domain  f2i  DD  fio  and  its  corresponding  domain 
of  dependence  Vfli . 

We  adopt  the  following  hypotheses. 

(i)  Approximation  properties  of  Mh  and  Pi,.  There  exist  constants  pm ,  sm,  with  0  <  [>m  <  f  and 
0  <  sm  <  r.  and  A m  such  that,  for  each  function  $  in  Cq°(Q o), 

\(uq  ~  PhUO,<p(0))\  <  Am  ||  U0  Hr.DQill  $  || Hi, 


li 


where  tp  is  the  solution  to  the  dual  problem  (3.12),  (3.13)  with  the  initial  data  $  for  the  dual  problem, 
(ii)  Residual.  Given  that  U  is  the  solution  to  (3.9),  (3.10),  there  exist  constants  pjv,  sjy.,  with  0  <  pjv  < 
£  and  0  <  sjv  <  r,  and  Ajv,  such  that  for  each  function  $  in  C'q0 (O0)  there  exists  x  G  G1  ([0, T]; A/*) 
with 


(Ut,tp  ~  x)h  +  B(U,  tp  -  \)  d t 


<  Aw  hp»+SN  ||  w0  ||  r,DQi  ||  ‘l5  || Hh 


where  tp  is  the  solution  to  the  dual  problem  (3.12),  (3.13)  with  the  initial  data  $  for  the  dual  problem, 
(iii)  Consistency.  Given  that  U  is  the  solution  to  (3.9),  (3.10),  there  exist  constants  sc  £  (0,  oo]  and 
A c  G  [0,  oo)  such  that 


B(U,tp)  d t 


<  A c  hSc  ||  Uq  Hr.'DOi 


$ 


U,Qo  j 


where  tp  is  the  solution  to  the  dual  problem  (3.12),  (3.13)  with  the  data  $  for  the  dual  problem. 
The  next  result  is  a  trivial  consequence  of  the  decomposition  (3.14)  and  conditions  (i)  -  (iii). 

Theorem  3.3.  Suppose  thatuo  G  L2per(Rd  ;Wn)CiHr  (Vfli;Wn) ,  with  Oo  CC  O i  C  I,  r  >  0,  and  assume 
that  conditions  (i)  -  (iii)  hold.  Then,  for  £  >  1,  we  have 


II  (u  —  U)(T)  ||  -ltQ0  <  C4  hs  ||  Uq  Jj^-DQi , 

where  s  =  min{pM  +  sm,  Pn  +  sjy,  }  and  C4  =  A m  +  A n  +  A c  • 

In  Table  3.2  we  display  the  parameters  of  the  above  result  for  some  finite  element  methods;  for  each  of 
the  methods  listed  we  have  Mh  =  A rh  and  have  assumed  that  Oi  =  I  (so  that  "DOi  =  I  also). 

Table  3.2 

The  parameters  of  Theorem  3.3  for  finite  element  methods  using  piecewise  polynomials  of  degree  k. 


parameter 

SG 

SG  with  AD 

DG 

Pm 

min{/c  +  1,  T} 

min-jA  +1,!} 

min{fc  +  1,  £} 

sm 

min{A:  +  1,  r} 

min{/c  +  1,  r} 

min{A;  +  1,  r} 

Pn 

min{/c,  £} 

min{A:,  /  } 

min{A;  +  1/2,  £} 

Sn 

min{A:,  r} 

min{A:,  r} 

inin{A  +  1/2,  r} 

Sc 

00 

7 

00 

3.3.  The  error  estimates.  Now  we  combine  the  results  obtained  in  the  previous  subsections. 
Theorem  3.4.  Letu  be  the  exact  solution  of  problem,  (3.4),  (3.5);  and  let  U  be  the  approximation  defined 
by  (3.9),  (3.10)  for  which  conditions  (i)  -  (iii)  are  valid.  Consider  the  convolution  kernel  K1’/1  of  Theorem 
3.1.  Let  each  of  the  components  of  u(T)  be  in  Hv(Tli)  and  let  Oo  be  such  that  Oo  +  2supp(A'G1)  CC  Oi. 
Then,  for  general  regular  triangulations  and  H  =  /js/T+fl  <  Hq,  we  have 

II  U(T)  —  K1^1  *  U(T)  |]o,o0  <Ch0s, 

where  0,  s  and  C  are  as  in  Theorem  3.3  with  C3  =  C4  |j  uo||r,'Pn1  and  6  =  v/(v  +  £). 

Moreover,  if  the  triangulation  is  translation  invariant  on  a  neighborhood  of  the  support  of  the  solution 
of  the  adjoint  equation  (3.12),  (3.13)  then,  for  H  =  h, 

II  u(T)  —  KVjjn  -k  U{T)  ||0,n0  <Chs, 
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C3  =  C4  ||  UoWr+t.vn-!  ■ 

Proof.  The  first  inequality  is  a  direct  consequence  of  Corollary  3.2  and  Theorem  3.3.  The  second 
inequality  also  follows  from  the  above  results  and  from  the  fact  that  if  the  triangulation  is  translation 
invariant  in  a  neighborhood  (of  order  H  =  h)  of  the  support  of  the  solution  of  the  adjoint  problem,  then  we 
have 


||  d°j{u  —  U)(T)  ||  -qo0  <  C4  hs  ||  d^uo  Ur/DOj- 

This  completes  the  proof.  □ 

Some  important  particular  cases  for  which  fl1  =  I  (and  consequently  Vfli  =  I)  are  collected  in  the  table 
below;  these  are  in  fact  the  estimates  we  can  actually  prove.  The  case  in  which  fh  ^  I  remains  a  challenging 
open  problem. 

Table  3.3 

Orders  of  convergence  with  piecewise  polynomials  of  degree  k  when  the  analytical  solution  u  is  in  C'([0,T];  H"er{I)). 


triangulations 

V 

SG 

SG  with  AD 

DG 

general 

0 

6  k 

6  a 

6(k  +  1/2) 

general 

2k  +  2 

9k 

9  a 

6(k  +  1/2) 

locally  invariant 

0 

k 

a 

k  +  1/2 

locally  invariant 

2A-  +  2 

2k 

a 

2k+l 

4.  Proofs. 

4.1.  The  approximation  result.  In  this  subsection,  for  the  sake  of  completeness,  we  sketch  the  proof 
of  Theorem  3.1  following  Bramble  and  Schatz  [4], 

Consider  the  following  quantity: 

Qh  =  II  u  -  Kff1  *  U  ||o,Q0  <  Qh,  1  +  Qh, 2, 

where 

Qh,i  —  \\u  —  i\h  *u||o,q0j 
®H,2  =  II  KV}fn  *  (u  —  U )  IIo.Qq  ■ 

To  estimate  Qh.  1,  we  denote  the  support  of  by  I,  we  label  the  Taylor  polynomial  of  degree  v  —  1  of 

u  around  y  by  Tvu{y ,  •),  and  we  put  1Zvu{y ,  •)  =  u(-)  —  Tvu{y,  •)•  We  then  easily  deduce  that 

u(x)  —  Kff1  *u(x)  =  Rvu(y,x)  —  J  I\I/'n(z )  Rvu(y,x  —  Hz)  dz, 

by  using  the  fact  that  the  kernel  K'ff1  reproduces  polynomials  of  degree  v  —  1  by  convolution.  For  y  =  x, 
the  above  expression  becomes 

u(x)  —  Kff1  ~k  u(x)  =  —  J  Kly,n(z)  Rvu{x,x  —  Hz)  dz, 

and  we  obtain 

Qh,  1  <  ||  Kv’n  ||Li(Rq  sup  ||  Rvu(-,  --Hz)  ||0,n0 

z£T 

Hv 

<  —  II  KV’1  Ili1  (X )  I  U  \v,Ua  +  HI- 
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On  applying  the  triangle  inequality  to  the  expression  of  Kv’n  given  by  (3.1),  we  get 

\\Kv’n\\LHX)<  E  I  K’n  I II  '0n_Q('  -  7)  IliPfR-q  =  E  \K’n\  =  c^ 

~,£Zd  ~,'£Zd 

since  ||  ipn~a(-  —  7)  ||Li(^)  =  1.  This  implies  that 

H1' 

Off,  1  <  Cl  I  u  |i,,Q0+ff  I  <  -j-y-  Cl  |  u  I^Qj . 

Now,  let  us  sketch  the  procedure  to  estimate  O h,i-  Take  a  set  2  such  that,  for  all  H  <  H0, 

O0  +  supp(A^)  C  fii/o, 
fii/2  +  supp(A^)  C  Sli . 


Then,  setting  e(x)  =  u(x )  —  U(x),  we  get 


0 


ff,  2 


|0,Qo  ^ 


<  C- 


Ei 

H<r 


D( 


'(K^n  *e)  \\->e,aiJ2, 


where  Co  depends  solely  on  f 20,  fli/o,  d,  j/,  and  f,  by  Lemma  4.2  in  Bramble  and  Schatz  [4].  This  is  the 
significant  step  that  allows  us  to  pass  from  the  A2-norm  to  a  negative-order  Sobolev  norm. 

Next,  we  exploit  the  fact  that  the  kernel  K^n  is  a  linear  combination  of  B-splines  given  by  (3.1);  this 
is  the  only  place  in  this  proof  where  properties  of  B-splines  are  used.  Thus,  by  the  property  (3.2)  we  have 
that 


D*(K%tl*e)  =  K%tl'’a*df[e, 

where 


Kv’tv‘«(y)  =  E  K’n^n-aHy-^). 

~,€Zd 


This  implies  that 


Off, 2  <  Co  E  \\K^m*d%eU,n1/2  <  C2  E  II  Kg1*  lUhK^II  \ 

\a\<(  \a\<t 

Finally,  since  ||  A)'/I:a  ||li(r4j  =  ||  ||Li^  <  C, .  we  get 

Off, 2  <  Ci  Co  E  II  <9ffe  ||-t,niy 
\m\<t 

This  completes  the  proof  of  Theorem  3.1. 


-e,Qi 


4.2.  The  conditions  (i)— (in)  for  some  finite  element  methods.  In  this  subsection,  we  justify  the 
results  displayed  in  Table  3.2. 

a.  The  SG  method.  Let  us  begin  by  considering  property  (i).  For  the  Zr-projection,  it  is  well  known 
that  pm  =  min{fc  +  1, £}  and  that  sm  =  min{fc  +  l,r}  for  regular  triangulations.  Next,  let  us  consider 
property  (iii).  A  simple  calculation  shows  that  we  can  take  AG  =  0;  this  allows  us  to  take  sc  =  00. 

For  property  (ii),  we  proceed  as  follows: 

0  =  /  {(Uti(p  ~  x)h  +  B(U,  ip  -  x)}  d t 
Jo 

{( ( U  -  u)t  +  A(U  -  u),tp  -  x)h}  d t. 
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Taking  y  as  the  Zr-projection  of  tp,  we  get  (ii)  with  p jv  =  min {k,(},  sn  =  min{fc,r}  and  =  I. 

b.  The  SG  with  artificial  diffusion.  For  this  method,  we  only  have  to  focus  on  property  (iii).  We 
have, 

m  <Pt)h  -  B(U,  p)}  d t  =  f  h?(VU,  Vp)  At 

Jo 

<  h 1  II  VG  Wl2(0,T-L2(I))\\  Vtp  Wl2(0,T;L2(I))- 


This  means  that  property  (iii)  is  satisfied  with  sc  =  7. 

c.  The  DG  method.  For  properties  (i)  and  (iii),  we  proceed  as  in  the  two  methods  discussed  above. 
The  verification  of  property  (ii)  requires  a  more  delicate  argument.  For  a  function  W  whose  components 
are  in  H^)erh{Wl\ Mm),  we  set  |.4W](a;)  =  A  ■  n+  W+{x)  +  A  ■  n~  W~{x)  for  every  x  G  e,  where  W^{x)  = 
lim~4_0  W(x  —  z  n ±)  and  uA  is  orthogonal  to  the  face  e  of  the  element  K  at  x.  With  this  notation,  we  can 
write  that 

0  =  /  {(Ut » V  ~  X)h  +  B(U,  ip  -  x)}  At 
Jo 

{Ut  +  AU.  p  -  x)h  +  [ A  Uj,tp  -  x)e  \  At, 

e€£  J 

where  we  obtained  the  last  step  after  a  simple  integration  by  parts;  by  £  we  denote  the  collection  of  faces  e 
of  the  elements  K  of  the  triangulation  Th.i- 

Now,  taking  y  to  be  the  T2-projection  of  tp  onto  Mh  =  A4,  we  get, 


This  implies  that  property  (ii)  is  satisfied  with  sn  =  min{fc  +  1/2,  £}  for  regular  triangulations.  It  remains 
to  obtain  an  estimate  of  the  first  term  on  the  right;  following  Johnson  and  Pitkaranta  [14]  or  Cockburn  [6], 
it  is  easy  to  prove  that  (iii)  is  satisfied  with  p n  =  niinjA:  +  1/2,  r}  and  =  I. 

5.  Numerical  experiments.  In  this  section,  we  validate  our  theoretical  results  with  an  emphasis  on 
the  case  in  which  the  doubling  of  the  order  of  convergence  is  achieved.  We  also  explore  the  performance 
of  the  post-processing  technique  in  situations  not  predicted  by  our  analysis;  thus,  we  display  the  1°° -errors 
in  all  our  experiments,  including  an  example  of  a  linear  convection-diffusion  equation  and  an  example  of  a 
non-linear  convection  equation. 

We  consider  the  discontinuous  Galerkin  method  with  polynomials  of  degree  k  (denoted  by  Pk)  and  use 
a  third  order  Runge-Kutta  method  to  discretize  in  time;  the  time  step  At  is  chosen  small  enough  so  that 
spatial  errors  dominate.  Results  for  P 1  to  P 4  are  shown. 

The  L°°  error  measures  the  maximum  numerical  error  at  the  six  Gaussian  points  in  each  element  for  all 
elements.  The  L2  error  is  computed  by  the  six-point  Gaussian  rule  in  each  element. 


Example  5.1.  A  linear  scalar  convection  equation  with  smooth  solution  on  the  domain  I  =  [0, 27r) : 

ut  +  ux  =0,  in  J  x  (0,T);  u(x,  0)  =  sin(®)  x  G  I, 


(5.1) 
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with  periodic  boundary  conditions.  The  errors  are  computed  at  T  =  12.5  which  is  about  2  periods  in  time. 

In  Table  5.1,  we  show  the  numerical  errors  for  this  problem.  We  can  clearly  see  that  both  the  L 2  and 
P°°  error  for  Pk  elements  is  of  ( k  +  l)-th  order  before  post-processing  and  of  at  least  (2k  +  l)-th  order  after 
post-processing.  This  is  consistent  with  our  theoretical  results. 


Table  5.1 

Example  5.1,  ut  +  ux  =  0,  smooth  solution. 


Before  postprocessing 

After  postprocessing 

mesh 

L2  error 

order 

P°°  error 

order 

L2  error 

order 

P°°  error 

order 

P1 

10 

3.29E-02 

— 

5.81E-02 

— 

3.01E-02 

— 

4.22E-02 

— 

20 

5.63E-03 

2.55 

1.06E-02 

2.45 

3.84E-03 

2.97 

5.44E-03 

2.96 

40 

1.16E-03 

2.28 

2.89E-03 

1.88 

4.79E-04 

3.00 

6.78E-04 

3.01 

80 

2.72E-04 

2.09 

8.08E-04 

1.84 

5.97E-05 

3.00 

8.45E-05 

3.00 

160 

6.68E-05 

2.03 

2.13E-04 

1.93 

7.45E-06 

3.00 

1.05E-05 

3.00 

320 

1.66E-05 

2.01 

5.45E-05 

1.96 

9.30E-07 

3.00 

1.32E-06 

3.00 

p2 

10 

8.63E-04 

— 

2.86E-03 

— 

2.52E-04 

— 

3.57E-04 

— 

20 

1.07E-04 

3.01 

3.69E-04 

2.95 

5.96E-06 

5.40 

8.41E-06 

5.41 

40 

1.34E-05 

3.00 

4.63E-05 

3.00 

1.53E-07 

5.29 

2.16E-07 

5.28 

80 

1.67E-06 

3.00 

5.78E-06 

3.00 

4.22E-09 

5.18 

5.97E-09 

5.18 

160 

2.09E-07 

3.00 

7.23E-07 

3.00 

1.27E-10 

5.06 

1.80E-10 

5.06 

p3 

10 

3.30E-05 

— 

9.59E-05 

— 

1.64E-05 

— 

2.31E-05 

— 

20 

2.06E-06 

4.00 

6.07E-06 

3.98 

7.07E-08 

7.85 

1.00E-07 

7.85 

40 

1.29E-07 

4.00 

3.80E-07 

4.00 

2.91E-10 

7.92 

4.15E-10 

7.91 

50 

5.29E-08 

4.00 

1.56E-07 

4.00 

5.03E-11 

7.87 

7.24E-11 

7.83 

p4 

10 

1.02E-06 

— 

2.30E-06 

— 

1.98E-06 

— 

2.81E-06 

— 

20 

3.21E-08 

5.00 

7.30E-08 

4.98 

2.20E-09 

9.82 

3.11E-09 

9.82 

30 

4.23E-09 

5.00 

9.66E-09 

4.99 

4.34E-11 

9.68 

6.66E-11 

9.48 

In  Fig.  5.1,  we  plot  the  errors  of  the  numerical  solution  before  and  after  post-processing  for  P 2  and  P3 
with  20  elements.  We  can  clearly  see  that  the  errors  before  post-processing  are  highly  oscillatory,  and  the 
post-processing  gets  rid  of  the  oscillation  in  the  error  and  greatly  reduces  its  magnitude. 

In  Fig.  5.2,  we  plot  the  errors,  in  absolute  value  and  in  logarithmic  scale,  of  the  numerical  solution 
before  and  after  post-processing  for  P2,  with  10,  20,  40,  80,  and  160  elements.  We  can  clearly  see  that 
the  post-processed  errors  are  less  oscillatory  and  much  smaller  in  magnitude,  and  approximately  third  and 
fifth  order  accuracy  for  the  pre-processed  and  post-processed  errors,  respectively,  measured  by  the  spacing 
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between  the  errors  when  the  number  of  elements  doubles. 


P2  with  20  elements,  Error 

before  post-processing  (solid  line)  ;  and  after  post-processing  (dashed  line) 


P3  with  20  elements,  Error 

before  post-processing  (solid  line)  ;  and  after  post-processing  (dashed  line) 


Fig.  5.1.  The  errors  before  and  after  post-processing  for  20  elements:  P 2  (left)  and  P3  (right). 


P2,  before  post-processing 


Fig.  5.2.  The  errors  in  absolute  value  and  in  logarithmic  scale,  for  P2  with  N=10,  20,  f0,  80,  and  160  elements.  Before 
post-processing  (left)  and  after  post-processing  (right). 

Example  5.2.  A  linear  scalar  convection  diffusion  equation  with  smooth  solution  on  the  domain  I  =  [0, 2i r): 

ut  +  ux  =  uxx,  in  /  x  (0 ,  T)  u(x,  0)  =  sin(x),  x  £  I,  (5.2) 

with  periodic  boundary  conditions.  The  errors  are  computed  at  T  =  2,  using  the  local  discontinuous  Galerkin 
method  [8],  with  alternating  left  and  right  fluxes  for  u  and  the  auxiliary  variable  q  which  approximates  ux 
(formula  (2.9)  in  [8]). 

Although  not  proven  in  this  paper,  we  expect  the  same  accuracy  result  to  hold  for  the  post-processed 
solution  as  in  Example  5.1.  In  Table  5.2,  we  show  the  results  for  this  problem.  We  can  clearly  see  that  the 
L°°  errors  for  Pk  elements  are  of  (k  +  l)-th  order  before  post-processing  and  of  at  least  (2k  +  l)-th  order 
after  post-processing,  both  for  the  solution  u  and  for  the  auxiliary  variable  q  which  approximates  ux.  The 
results  for  the  L 2  errors  are  similar  and  are  not  shown  to  save  space. 

Example  5.3.  The  same  linear  scalar  convection  equation  (5.1)  with  the  same  initial  condition,  except  that 
now  I  =  [0,5).  The  solution  now  has  a  discontinuity  at  x  =  0  (or  x  =  5)  and  this  discontinuity  moves  in 
time  with  the  characteristic  speed  1.  We  compute  the  errors  at  T  =  12.5,  that  is,  after  2.5  periods  in  time. 
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Table  5.2 

Example  5.2,  ut  +  ux  =  uxx,  smooth  solution. 


Before  postprocessing 

After  postprocessing 

u 

q  for  ux 

u 

q  for  ux 

mesh 

L°°  error 

order 

poc  error 

order 

poc  error 

order 

poo  error 

order 

P 1 

10 

6.74E-03 

— 

5.82E-03 

— 

1.19E-03 

— 

1.18E-03 

— 

20 

1.82E-03 

1.89 

1.71E-03 

1.77 

1.34E-04 

3.15 

1.41E-04 

3.07 

40 

4.68E-04 

1.96 

4.56E-04 

1.91 

1.56E-05 

3.05 

1.69E-05 

3.06 

80 

1.19E-04 

1.98 

1.17E-04 

1.96 

1.46E-06 

3.02 

2.07E-06 

3.03 

160 

3.00E-05 

1.99 

2.98E-05 

1.98 

2.32E-07 

3.03 

2.57E-07 

3.01 

320 

7.52E-06 

1.99 

7.50E-06 

1.99 

2.87E-08 

3.01 

3.20E-08 

3.01 

P 2 

10 

3.97E-04 

— 

3.38E-04 

— 

2.93E-05 

— 

2.96E-05 

— 

20 

5.01E-05 

2.99 

4.61E-05 

2.87 

5.43E-07 

5.75 

5.46E-07 

5.76 

40 

6.25E-06 

3.00 

6.02E-06 

2.94 

1.04E-08 

5.71 

1.05E-08 

5.70 

80 

7.83E-07 

3.00 

7.68E-07 

2.97 

2.19E-10 

5.57 

2.26E-10 

5.54 

160 

9.78E-08 

3.00 

9.69E-08 

2.99 

5.31E-12 

5.37 

5.63E-12 

5.33 

p3 

10 

1.30E-05 

— 

1.13E-05 

— 

3.09E-06 

— 

3.09E-06 

— 

20 

8.23E-07 

3.98 

7.73E-07 

3.86 

1.32E-08 

7.87 

1.32E-08 

7.87 

40 

5.14E-08 

4.00 

4.99E-08 

3.95 

5.31E-11 

7.96 

5.31E-11 

7.96 

pA 

10 

3.11E-07 

— 

2.93E-07 

— 

3.79E-07 

— 

3.80E-07 

— 

20 

9.89E-09 

4.97 

9.54E-09 

4.94 

4.19E-10 

9.82 

4.19E-10 

9.82 

30 

1.30E-09 

5.00 

1.27E-09 

4.98 

7.11E-12 

10.05 

7.11E-12 

10.05 

The  discontinuity  at  this  time  is  located  at  x  =  2.5.  The  errors  shown  in  Table  5.3  are  calculated  within 
the  “smooth  region”  [0, 1]  U  [4,5],  at  distance  1.5  away  from  the  discontinuity,  namely  excluding  the  interval 
1  <  x  <  4. 

The  theory  in  this  paper  would  only  guarantee  ( k  +  l)-th  order  accuracy  for  Pk  elements  after  post¬ 
processing  since  our  estimates  hold  for  Vfli  =  I  only  and  the  initial  condition  displays  a  discontinuity. 
However,  Table  5.3  shows  that  both  the  L 2  errors  and  the  L°°  errors  are  still  at  least  (2k  +  l)-th  order 
accurate  for  Pk  elements  after  post-processing.  This  indicates  that  it  is  reasonable  to  expect  that  a  similar 
result  with  a  domain  DOi  excluding  the  discontinuity  should  hold,  see  Fig.  5.4. 

In  Fig.  5.3  we  plot  the  errors,  in  absolute  value  and  in  logarithmic  scale,  of  the  numerical  solution 
before  and  after  post-processing  for  P 2,  with  10,  20,  40,  80  and  160  elements.  We  can  clearly  see  that  the 
post-processed  errors  are  less  oscillatory  and  much  smaller  in  magnitude  away  from  the  discontinuity. 
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Table  5.3 

Example  5.3,  ut  +  ux  =  0,  discontinuous  solution,  errors  in  smooth  regions 


Before  postprocessing 

After  postprocessing 

mesh 

L2  error 

order 

P°°  error 

order 

L2  error 

order 

P°°  error 

order 

P1 

10 

2.02E-02 

— 

6.46E-02 

— 

1.76E-02 

— 

2.80E-02 

— 

20 

4.37E-03 

2.21 

1.21E-02 

2.41 

3.96E-03 

2.15 

1.18E-02 

1.24 

40 

6.63E-04 

2.72 

1.89E-03 

2.69 

2.69E-04 

3.88 

6.77E-04 

4.12 

80 

1.58E-04 

2.07 

5.24E-04 

1.85 

2.78E-05 

3.27 

4.31E-05 

3.97 

160 

3.92E-05 

2.01 

1.36E-04 

1.94 

3.49E-06 

3.00 

5.31E-06 

3.02 

320 

9.80E-06 

2.00 

3.47E-05 

1.97 

4.37E-07 

3.00 

6.63E-07 

3.00 

p2 

10 

4.53E-03 

— 

1.08E-02 

— 

3.74E-03 

— 

1.15E-02 

— 

20 

4.96E-04 

3.19 

1.98E-03 

2.45 

3.02E-04 

3.63 

1.07E-03 

3.42 

40 

8.80E-06 

5.82 

2.51E-05 

6.30 

4.03E-06 

6.23 

2.74E-05 

5.29 

80 

8.97E-07 

3.29 

2.91E-06 

3.11 

1.74E-09 

11.18 

1.32E-08 

11.02 

160 

1.12E-07 

3.00 

3.64E-07 

3.00 

5.09E-11 

5.09 

8.75E-11 

7.23 

p3 

10 

2.87E-03 

— 

1.24E-02 

— 

7.76E-04 

— 

1.81E-03 

— 

20 

1.97E-04 

3.87 

1.03E-03 

3.60 

6.91E-06 

6.81 

2.92E-05 

5.95 

40 

1.36E-06 

7.18 

7.21E-06 

7.15 

3.51E-08 

7.62 

1.88E-07 

7.27 

80 

3.03E-09 

8.81 

1.01E-08 

9.47 

2.18E-11 

10.65 

6.89E-11 

11.42 

p4 

10 

1.93E-03 

— 

6.32E-03 

— 

1.36E-03 

— 

2.91E-03 

— 

20 

9.79E-05 

4.30 

5.42E-04 

3.54 

1.15E-07 

13.53 

8.37E-07 

11.76 

40 

5.86E-07 

7.39 

4.70E-06 

6.85 

3.46E-11 

11.70 

2.11E-10 

11.96 

Example  5.4.  A  scalar  nonlinear  Burgers’  equation  with  continuous  and  discontinuous  solutions  on  the 
domain  I  =  [0, 2ir): 


ut  + 


in  I  x  (0, T); 


u(x,  0)  =  -  +  sin(x),  x  G  I] 


(5.3) 


with  periodic  boundary  conditions.  The  errors  at  T  =  0.5,  when  the  solution  is  still  smooth,  are  given  in 
Table  5.4.  It  seems  that  in  general,  post-processed  errors  are  still  smaller,  although  the  asymptotic  orders 
seem  to  show  up  later  than  for  the  linear  case,  as  the  mesh  is  refined.  We  remark  that  the  theory  in  this 
paper  does  not  cover  this  nonlinear  problem. 

In  Fig.  5.5,  we  plot  the  errors  of  the  numerical  solution  before  and  after  post-processing  for  P 2  and  P3 
with  20  elements.  From  Table  5.4  we  can  see  that  in  both  situations  the  errors  after  post-processing  are 
actually  larger  than  the  errors  before  post-processing.  Note  in  Fig.  5.5  that  near  the  middle  region,  the 
oscillations  in  the  errors  are  not  “uniform”,  apparently  due  to  nonlinear  effects,  hence  the  post-processing 
actually  gives  larger  errors.  Fortunately,  for  a  larger  number  of  elements  the  post-processing  begins  to  be 
effective  and  the  errors  after  post-processing  do  become  smaller,  see  Table  5.4  and  the  following  Fig.  5.6. 
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In  Fig.  5.6  we  plot  the  errors,  in  absolute  value  and  in  logarithmic  scale,  of  the  numerical  solution  before 
and  after  post-processing  for  P2 ,  with  10,  20,  40,  80,  160  and  320  elements.  We  can  clearly  see  that  the 
post-processed  errors  are  less  oscillatory  and  much  smaller  in  magnitude,  especially  for  a  large  number  of 
elements.  However,  notice  that  due  to  non-linear  effects  not  all  oscillations  in  the  errors  have  been  removed 
by  the  post-processing,  especially  for  a  large  number  of  elements. 

Table  5.4 

Example  5-4,  Burgers’  equation  with  smooth  solution. 


Before  postprocessing 

After  postprocessing 

mesh 

L 2  error 

order 

L°°  error 

order 

L 2  error 

order 

L°°  error 

order 

P 1 

10 

1.95E-02 

— 

8.87E-02 

— 

1.37E-02 

— 

3.99E-02 

— 

20 

5.31E-03 

1.88 

2.77E-02 

1.68 

1.63E-03 

3.07 

6.47E-03 

2.63 

40 

1.33E-03 

2.00 

7.55E-03 

1.87 

1.28E-04 

3.67 

5.55E-04 

3.54 

80 

3.33E-04 

1.99 

1.95E-03 

1.95 

1.03E-05 

3.63 

4.17E-05 

3.73 

160 

8.37E-05 

1.99 

4.99E-04 

1.97 

1.12E-06 

3.20 

4.21E-06 

3.31 

320 

2.10E-05 

2.00 

1.26E-04 

1.98 

1.42E-07 

2.98 

5.69E-07 

2.89 

P 2 

10 

3.46E-03 

— 

1.93E-02 

— 

1.12E-02 

— 

3.37E-02 

— 

20 

4.81E-04 

2.85 

3.57E  03 

2.43 

9.25E-04 

3.59 

3.47E-03 

3.28 

40 

8.00E-05 

2.59 

6.22E-04 

2.52 

3.63E-05 

4.67 

1.58E-04 

4.46 

80 

1.30E-05 

2.62 

1.20E-04 

2.37 

8.43E-07 

5.43 

3.93E-06 

5.33 

160 

2.04E-06 

2.67 

1.98E-05 

2.61 

1.67E-08 

5.66 

8.51E-08 

5.53 

320 

3.06E-07 

2.73 

3.02E-06 

2.71 

3.60E-10 

5.53 

1.85E-09 

5.52 

p3 

10 

4.33E-04 

— 

2.24E-03 

— 

1.12E-02 

— 

3.35E-02 

— 

20 

4.16E-05 

3.38 

2.00E-04 

3.48 

8.08E-04 

3.80 

3.03E-03 

3.46 

40 

2.43E-06 

4.10 

1.74E-05 

3.53 

2.06E-05 

5.30 

9.42E-05 

5.01 

80 

1.46E-07 

4.06 

1.04E-06 

4.07 

1.96E-07 

6.72 

1.01E-06 

6.54 

160 

1.03E-08 

3.82 

6.72E-08 

3.95 

1.10E-09 

7.47 

5.94E-09 

7.41 

P 4 

10 

1.75E-04 

— 

8.25E-04 

— 

1.15E-02 

— 

3.36E-02 

— 

20 

4.19E-06 

5.39 

2.45E-05 

5.07 

7.63E-04 

3.91 

2.82E-03 

3.58 

40 

1.70E-07 

4.62 

1.04E-06 

4.55 

1.48E-05 

5.69 

6.82E-05 

5.37 

50 

6.45E-08 

4.36 

4.40E-07 

3.87 

3.09E-06 

7.03 

1.52E-05 

6.74 

Next,  we  compute  the  solution  at  T  =  2,  that  is,  after  the  shock  has  developed.  We  measure  the  errors 
on  the  smooth  region  0.57T  away  from  the  discontinuity  and  show  the  results  in  Table  5.5.  The  codes  ran 
stably  only  for  P1  and  P2;  hence  only  these  two  cases  are  shown. 

In  order  to  stabilize  the  algorithm  in  the  presence  of  shocks,  we  apply  a  TVB  (total  variation  bounded) 
limiter  with  M  =  3,  see  [7].  This  limiter  has  no  effect  on  the  numerical  solution  at  T  =  0.5  when  the  solution 
is  still  smooth,  but  allows  the  scheme  to  run  stably  for  P3  and  P 4  after  the  shock  develops.  We  again 
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Table  5.5 

Example  5..(,  Burgers’  equation  with  discontinuous  solution. 


Before  postprocessing 

After  postprocessing 

mesh 

L2  error 

order 

pec  error 

order 

L 2  error 

order 

L°°  error 

order 

P 1 

10 

8.70E-03 

— 

3.56E-02 

— 

6.79E-03 

— 

1.99E-02 

— 

20 

3.05E-04 

4.83 

1.47E-03 

4.60 

2.23E-04 

4.93 

8.61E-04 

4.53 

40 

1.70E-05 

4.16 

8.14E-05 

4.18 

1.09E-05 

4.36 

2.25E-05 

5.26 

80 

3.71E-06 

2.20 

2.07E-05 

1.97 

1.37E-06 

2.99 

2.86E-06 

2.97 

160 

8.65E-07 

2.10 

4.67E-06 

2.15 

1.63E-07 

3.07 

3.43E-07 

3.06 

320 

2.17E-07 

2.00 

1.19E-06 

1.97 

2.05E-08 

3.00 

4.31E-08 

2.99 

P 2 

10 

6.26E-03 

— 

3.29E-02 

— 

1.57E-03 

— 

7.05E-03 

— 

20 

2.77E-04 

4.50 

1.44E-03 

4.52 

5.47E-05 

4.84 

1.52E-04 

5.54 

40 

2.03E-05 

3.77 

1.68E-04 

3.10 

6.88E-06 

2.99 

2.62E-05 

2.53 

80 

2.30E-06 

3.14 

2.17E-05 

2.95 

8.39E-07 

3.03 

4.61E-06 

2.51 

160 

4.23E-07 

2.44 

4.75E-06 

2.19 

1.16E-07 

2.86 

7.95E-07 

2.54 

320 

6.12E-08 

2.79 

7.77E-07 

2.61 

1.41E-08 

3.04 

1.29E-07 

2.62 

measure  the  errors  on  the  smooth  region  0.57 r  away  from  the  discontinuity  and  show  the  result  in  Table  5.6. 

In  Fig.  5.7  we  plot  the  errors,  in  absolute  value  and  in  logarithmic  scale,  of  the  numerical  solution  before 
and  after  post-processing  for  P 2  with  a  TVB  limiter,  at  t  =  2,  with  10,  20,  40,  80,  160  and  320  elements.  We 
can  clearly  see  that  the  post-processed  errors  are  less  oscillatory  and  much  smaller  in  magnitude,  especially 
for  a  large  number  of  elements,  away  from  the  discontinuity.  Again,  notice  that  not  all  oscillations  in  the 
errors  have  been  removed  by  the  post-processing,  especially  for  a  large  number  of  elements,  due  to  non-linear 
effects. 


Example  5.5.  A  linear  system  with  a  smooth  solution  in  the  domain  I  =  [0,  2tt) : 


ut  +  vx  =  0 
vt  +  ux  =  0 


in  I  x  (0,  T), 


u(x,  0)  =  sin  (a) 
v(x,  0)  =  0 


X  €  I, 


(5.4) 


with  periodic  boundary  conditions.  The  errors  are  computed  at  T  =  12.5  which  is  about  2  periods  in  time. 

In  Table  5.7,  we  show  the  results  for  this  problem.  The  errors  are  the  combined  ones  of  u  and  v.  We  can 
clearly  see  that  both  L 2  and  L°°  errors  for  Pk  elements  are  (k  +  l)-th  order  before  post-processing  and  at 
least  (2k  +  l)-th  order  after  post-processing.  In  fact,  the  errors  are  very  similar  to  the  scalar  case  in  Example 
5.1.  This  is  consistent  with  our  theoretical  results. 

Notice  that  this  example  and  the  next  one  with  discontinuous  solution  for  linear  systems  indicate  that 
the  method  is  very  suitable  for  long  time  simulation  of  linear  systems  as  the  post-processing  needs  to  be 
performed  only  at  the  final  time.  Examples  include  aeroacoustic  problems  when  linear  Euler  equations  must 
be  solved  for  a  long  time  to  propagate  the  pressure  waves. 


Example  5.6.  The  same  linear  system  (5.4)  with  the  same  initial  condition,  except  that  now  0  <  x  <  5 
and  the  boundary  condition  is  5-periodic.  The  solution  now  has  a  discontinuity  at  x  =  0  (or  x  =  5)  and  this 
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Table  5.6 

Example  5-4,  Burgers’  equation  with  discontinuous  solution.  TVB  limiters. 


Before  postprocessing 

After  postprocessing 

mesh 

L 2  error 

order 

poc  error 

order 

L 2  error 

order 

pec  error 

order 

P 1 

10 

1.26E-03 

— 

4.44E-03 

— 

1.02E-03 

— 

2.28E-03 

— 

20 

1.16E-04 

3.44 

4.38E-04 

3.34 

1.01E-04 

3.33 

1.94E-04 

3.55 

40 

1.72E-05 

2.76 

8.33E-05 

2.40 

1.09E-05 

3.22 

2.26E-05 

3.11 

80 

3.72E-06 

2.20 

2.08E-05 

2.00 

1.37E-06 

2.99 

2.87E-06 

2.98 

160 

8.73E-07 

2.09 

4.75E-06 

2.13 

1.63E-07 

3.07 

3.44E-07 

3.06 

320 

2.17E-07 

2.01 

1.19E-06 

2.00 

2.05E-08 

3.00 

4.32E-08 

2.99 

P 2 

10 

1.03E-02 

— 

5.83E-02 

— 

1.99E-03 

— 

7.66E-03 

— 

20 

4.22E-04 

4.60 

3.36E-03 

4.12 

5.32E-05 

5.22 

1.50E-04 

5.68 

40 

2.23E-05 

4.24 

1.99E-04 

4.08 

6.87E-06 

2.95 

2.53E-05 

2.53 

80 

2.29E-06 

3.28 

2.16E-05 

3.20 

8.39E-07 

3.03 

4.61E-06 

2.50 

160 

4.21E-07 

2.44 

4.72E-06 

2.19 

1.16E-07 

2.86 

7.95E-07 

2.54 

320 

6.10E-08 

2.79 

7.75E-07 

2.61 

1.41E-08 

3.04 

1.29E-07 

2.62 

p3 

10 

9.98E-04 

— 

6.64E-03 

— 

3.45E-03 

— 

1.35E-02 

— 

20 

1.47E-04 

2.76 

1.38E-03 

2.20 

9.35E-06 

8.52 

5.18E-05 

8.03 

40 

4.92E-07 

8.22 

5.43E-06 

7.99 

2.92E-08 

8.32 

2.08E-07 

7.96 

80 

8.58E-10 

9.16 

1.43E-08 

8.57 

3.7  IE  10 

6.30 

8.73E-10 

7.90 

pA 

10 

3.86E-01 

— 

9.90E-01 

— 

2.28E-01 

— 

3.78E-01 

— 

20 

1.08E-01 

1.84 

2.16E-01 

2.20 

5.28E-02 

2.11 

1.34E-01 

1.50 

40 

1.89E-03 

5.83 

2.25E-02 

3.26 

3.88E-04 

7.09 

3.97E-03 

5.07 

80 

8.08E-08 

14.52 

6.11E-07 

15.17 

1.46E-08 

14.70 

7.42E-08 

15.71 

discontinuity  moves  in  time  with  the  characteristic  speeds  ±1.  We  compute  the  errors  at  t  =  12.5,  after  2.5 
periods  in  time.  The  two  discontinuities  at  this  time  are  both  located  at  x  =  2.5.  The  errors  shown  in  Table 
5.8  are  calculated  within  the  “smooth  region”  that  lies  a  distance  1.5  away  from  the  discontinuities,  namely 
excluding  the  interval  1  <  x  <  4. 

The  theory  in  this  paper  would  only  guarantee  ( k  +  l)-th  order  accuracy  for  Pk  elements  after  post¬ 
processing,  since  when  we  take  VQ i  =  I,  the  initial  condition  has  a  discontinuity  in  this  set.  However,  Table 
5.8  shows  that  both  the  L'2  errors  and  the  L°°  errors  are  still  at  least  {2k  +  l)-th  order  accurate  for  Pk 
elements  after  post-processing. 

In  fact,  the  behavior  of  the  errors  is  very  similar  to  that  observed  in  the  scalar  case  in  Example  5.3.  This 
is  not  really  surprising  since  our  linear  system  is  equivalent  to  the  following  two  decoupled  scalar  equations: 


( u  -  v)t  -  (u  -  v)x  =  0 
(u  +  v)t  +  { u  +  v)x  =  0 


in  /  x  (0, T), 
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( u  —  v){x ,  0)  =  sin(;r) 
( u  +  v)(x,  0)  =  sin(;c) 


x  €  I. 


Table  5.7 

Example  5.5,  linear  system  with  smooth  solution 


Before  postprocessing 

After  postprocessing 

mesh 

L2  error 

order 

pec  error 

order 

L 2  error 

order 

pec  error 

order 

P 1 

10 

2.33E-02 

— 

5.20E-02 

— 

2.13E-02 

— 

4.16E-02 

— 

20 

3.98E-03 

2.55 

8.55E-03 

2.60 

2.72E-03 

2.97 

5.36E-03 

2.96 

40 

8.20E-04 

2.28 

1.89E-03 

2.18 

3.39E-04 

3.00 

6.74E-04 

2.99 

80 

1.92E-04 

2.09 

4.77E-04 

1.98 

4.23E-05 

3.00 

8.43E-05 

3.00 

160 

4.72E-05 

2.03 

1.20E-04 

2.00 

5.28E-06 

3.00 

1.05E-05 

3.00 

320 

1.17E-05 

2.01 

2.99E-05 

2.00 

6.59E-07 

3.00 

1.31  E  06 

3.00 

P 2 

10 

6.10E-04 

— 

1.67E-03 

— 

1.78E-04 

— 

3.53E-04 

— 

20 

7.57E-05 

3.01 

2.08E-04 

3.00 

4.22E-06 

5.40 

8.42E-06 

5.39 

40 

9.46E-06 

3.00 

2.60E-05 

3.00 

1.09E-07 

5.28 

2.17E-07 

5.28 

80 

1.18E-06 

3.00 

3.24E-06 

3.00 

3.11E-09 

5.13 

6.20E-09 

5.13 

160 

1.48E-07 

3.00 

4.06E-07 

3.00 

8.95E-11 

5.12 

1.77E-10 

5.13 

p3 

10 

2.33E-05 

— 

5.73E-05 

— 

1.16E-05 

— 

2.30E-05 

— 

20 

1.46E-06 

4.00 

3.61E-06 

3.99 

5.00E-08 

7.85 

9.98E-08 

7.85 

40 

9.13E-08 

4.00 

2.27E-07 

3.99 

2.13E-10 

7.88 

4.25E-10 

7.88 

50 

3.74E-08 

4.00 

9.29E-08 

4.01 

3.94E-11 

7.56 

7.84E-11 

7.57 

P4 

10 

7.24E-07 

— 

1.37E-06 

— 

1.40E-06 

— 

2.79E-06 

— 

20 

2.27E-08 

5.00 

4.33E-08 

4.99 

1.56E-09 

9.82 

3.11E-09 

9.81 

30 

2.99E-09 

5.00 

5.72E-09 

4.99 

3.06E-11 

9.69 

6.16E-11 

9.67 

As  a  consequence,  the  domain  of  dependence  T>Qi  does  not  include  the  discontinuity  of  the  initial  condition; 
see  Fig.  5.8  (top). 

On  the  other  hand,  it  is  interesting  to  point  out  that  this  doubling  of  the  order  of  accuracy  does  not 
take  place  for  the  problem, 


Uf  —  ux  =  v 
vt  +  vx  =  -u 


in  I  x  (0 ,T), 


u(x,  0)  =  sin(x) 
v(x,  0)  =  0 


x  E  /, 


(5.5) 


with  periodic  boundary  conditions,  since  now  the  two  scalar  equations  associated  with  the  diagonalization  of 
the  system  are  coupled  through  zero-order  terms;  as  a  consequence,  the  domain  of  dependence  Df h  always 
includes  the  discontinuity  of  the  initial  condition;  see  Fig  5.8  (bottom).  This  is  the  example  treated  in 
the  early  work  of  Majda  and  Osher  and  [18]  and  Majda,  McDonough  and  Osher  [17].  Due  to  this  lack  of 
regularity  of  the  initial  condition  on  Vfli ,  post-processing  with  a  kernel  of  support  of  order  h  does  not  yield 
any  significant  improvement;  a  kernel  of  support  of  order  almost  one  is  required,  as  predicted  by  our  main 
theorem;  see  also  Mock  and  Lax  [19]. 
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Table  5.S 

Example  5.6,  linear  system  with  discontinuous  solution 


Before  postprocessing 

After  postprocessing 

mesh 

L2  error 

order 

pec  error 

order 

L 2  error 

order 

pec  error 

order 

P 1 

10 

1.49E-02 

— 

4.00E-02 

— 

1.29E-02 

— 

4.27E-02 

— 

20 

3.19E-03 

2.22 

9.35E-03 

2.10 

2.79E-03 

2.21 

7.69E-03 

2.47 

40 

4.76E-04 

2.74 

1.37E-03 

2.78 

1.85E-04 

3.91 

4.82E-04 

3.99 

80 

1.13E-04 

2.08 

3.04E-04 

2.17 

1.99E-05 

3.22 

4.28E-05 

3.49 

160 

2.78E-05 

2.02 

7.59E-05 

2.00 

2.48E-06 

3.00 

5.31E-06 

3.01 

320 

6.94E-06 

2.00 

1.90E-05 

2.00 

3.09E-07 

3.00 

6.63E-07 

3.00 

P 2 

10 

3.41E-03 

— 

6.93E-03 

— 

2.65E-03 

— 

9.29E-03 

— 

20 

3.58E-04 

3.25 

1.30E-03 

2.42 

2.22E-04 

3.58 

6.40E-04 

3.86 

40 

6.30E-06 

5.83 

2.42E-05 

5.74 

2.85E-06 

6.28 

1.39E-05 

5.52 

80 

6.33E-07 

3.32 

1.64E-06 

3.88 

1.23E-09 

11.17 

7.52E-09 

10.86 

160 

7.91E-08 

3.00 

2.05E-07 

3.00 

3.34E-11 

5.20 

5.54E-11 

7.08 

p3 

10 

2.03E-03 

— 

6.43E-03 

— 

5.35E-04 

— 

1.77E-03 

— 

20 

1.40E-04 

3.86 

5.41E-04 

3.57 

4.92E-06 

6.76 

2.27E-05 

6.28 

40 

9.66E-07 

7.18 

3.64E-06 

7.21 

2.50E-08 

7.62 

9.60E-08 

7.89 

80 

2.14E-09 

8.82 

6.00E-09 

9.25 

1 .34  E  1 1 

10.87 

4.89E-11 

10.94 

pA 

10 

1.38E-03 

— 

3.26E-03 

— 

9.61E-04 

— 

2.90E-03 

— 

20 

6.92E-05 

4.32 

2.72E-04 

3.58 

8.09E-08 

13.54 

6.31E-07 

12.17 

40 

4.14E-07 

7.39 

2.36E-06 

6.85 

2.42E-11 

11.71 

1.03E-10 

12.58 

6.  Extensions  and  concluding  remarks.  In  this  paper,  we  have  shown  how  to  enhance  the  approxi¬ 
mation  given  by  a  finite  element  method  for  linear  hyperbolic  equations  by  applying  a  simple  post-processing 
at  the  very  end  of  the  computations.  Our  theoretical  results  can  be  easily  extended  to  the  case  in  which  the 
matrices  Aj,j  =  0 ,d,  are  very  smooth  functions  of  (x,t).  To  do  that,  it  is  enough  to  mimic  the  induction 
argument  presented  by  Bramble  and  Schatz  in  [4]. 

The  role  of  negative-order  error  estimates  is  crucial  since  it  is  the  analytical  tool  that  captures  the 
ocillatory  nature  of  the  error.  For  these  negative-order  norms  of  the  error,  upper  bounds  were  obtained 
which  depend  on  a  global  norm  of  the  initial  data.  Our  numerical  results  suggest,  however,  that  they  should 
depend  only  on  a  local  norm  of  the  initial  data.  In  fact,  a  result  of  this  type  was  obtained  in  1998  for  finite 
difference  schemes  by  Engquist  and  Sjogreen  [10].  To  obtain  such  a  result  for,  say,  the  discontinuous  Galerkin 
method  is  an  challenging  open  problem. 

Finally,  let  us  end  by  pointing  out  that  our  numerical  results  seem  to  indicate  that  the  post-processing 
has  a  positive  impact  on  the  quality  of  the  approximate  solution  even  if  the  problem  is  non-linear.  A 
theoretical  analysis  of  this  case  is  yet  another  important  open  problem. 
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Fig.  5.6.  The  errors  in  absolute  value  and  in  logarithmic  scale,  for  P2  with  N=20,  40,  80,  160  and  320  elements.  Smooth 
solution  of  Burgers  equation.  Before  post-processing  (left)  and  after  post-processing  (right). 


Fig.  5.7.  The  errors  in  absolute  value  and  in  logarithmic  scale,  for  P 2  with  N=10,  20,  )0,  80,  160  and  320  elements. 
Discontinuous  solution  of  Burgers  equation.  Before  post-processing  (left)  and  after  post-processing  (right). 
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Fig.  5.8.  The  domain  of  smoothness  of  u(T),  fig,  the  domain  f!i  DD  Qo  (and  its  corresponding  domain  of  dependence 
VTli  for  the  system  (5.4)  (top)  and  the  system  (5.5)  (bottom).  Note  the  discontinuity  curves  t=  \x\. 
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